/*    Copyright 2010 Tobias Marschall
 *
 *    This file is part of MoSDi.
 *
 *    MoSDi is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    MoSDi is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoSDi.  If not, see <http://www.gnu.org/licenses/>.
 */

package mosdi.subcommands;

import java.io.BufferedReader;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.IOException;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.StringTokenizer;

import mosdi.distributions.InhomogeneousDistribution;
import mosdi.util.Iupac;
import mosdi.util.Log;
import mosdi.util.LogSpace;

public class CondPvalueSubcommand extends Subcommand {
	
	@Override
	public String usage() {
		return
		super.usage()+" [options] <inhom-distribution-file> <inputfile>\n" +
		"\n" +
		"Options:\n" +
		"  -t <pvalue-threshold>: dont print pattern if conditional pvalue is above threshold\n" +
		"  -r: simultaneously consider reverse complementary motif";
	}
	
	@Override
	public String description() {
		return 
		"Takes the output of \""+(new DiscoverySubcommand().name())+"\" as input and computes the p-values " +
		"conditional on the expectation w.r.t. a given inhomogeneous model.";
	}

	@Override
	public String name() {
		return "cond-pvalue";
	}

	@Override
	public int run(String[] args) {
		parseOptions(args, 2, "t:r");

		// Option dependencies
		// -- none --

		// Mandatory arguments
		String inhomDistFilename = getStringArgument(0);
		String inputFilename = getStringArgument(1);

		// Options
		double pvalueThreshold = getRangedDoubleOption("t", 0.0, 1.0, 1.0);
		boolean considerReverse = getBooleanOption("r", false);
		
		InhomogeneousDistribution inhomDist = InhomogeneousDistribution.readFromFile(inhomDistFilename);
			
		InhomogeneousDistribution.QGram[] inhomDistQGrams = null;
		int inhomDistQGramsQ = -1;
		FileInputStream file = null;
		try {
			file = new FileInputStream(inputFilename);
		} catch (FileNotFoundException e) {
			Log.errorln("Inputfile not found, sorry!");
			System.exit(1);
		}
		BufferedReader br = new BufferedReader(new InputStreamReader(file));
		try {
			while (true) {
				String line = br.readLine();
				if (line==null) break;
				// tokenize ...
				StringTokenizer st = new StringTokenizer(line," ",false);
				ArrayList<String> tokens = new ArrayList<String>();
				while (st.hasMoreTokens()) tokens.add(st.nextToken());
				// check ...
				int matches = -1;
				try {
					if (tokens.size()<9) throw new IllegalArgumentException();
					matches = Integer.parseInt(tokens.get(8));
					if (tokens.size()!=19+matches) throw new IllegalArgumentException();
					if (!tokens.get(0).equals(">>p_value>>")) throw new IllegalArgumentException();
					if (!tokens.get(17).equals(">>p_value_table>>")) throw new IllegalArgumentException();
				} catch (IllegalArgumentException e) {
					Log.errorln("Invalid line:");
					Log.errorln(line);
					System.exit(1);
				}
				String pattern = tokens.get(4);
				
					
				if (inhomDistQGramsQ!=pattern.length()) {
					Log.startTimer();
					inhomDistQGrams = inhomDist.getQGrams(pattern.length());
					Log.stopTimer("Inhomogeneous model: create qgram map");
					inhomDistQGramsQ = pattern.length();
				}
				
				Log.startTimer();
				double expectation = 0.0;
				expectation+=inhomDist.calcExpectation(Iupac.toGeneralizedString(pattern),inhomDistQGrams);
				if (considerReverse) {
					String reversePattern = Iupac.reverseComplementary(pattern);
					expectation+=inhomDist.calcExpectation(Iupac.toGeneralizedString(reversePattern),inhomDistQGrams);
				}
				Log.stopTimer("Calculation of expectation based on inhomogeneous model");
				Log.printf(Log.Level.VERBOSE, "Expectation (inhomogeneous null model): %e%n", expectation);
				if (expectation>=matches) {
					Log.println(Log.Level.STANDARD, "Occurrences below expectation --> skipping");
					continue;
				}
				boolean logarithmic = tokens.get(2).equals("LOG");
				double orginalPValue = Double.parseDouble(tokens.get(18 + matches));
				double expPValue = Double.parseDouble(tokens.get(18 + (int)Math.round(expectation)));
				if (logarithmic) {
					double newPValue = orginalPValue-expPValue;
					Log.printf(Log.Level.STANDARD, "former pvalue: %s, exp. pvalue: %s --> pvalue: %s%n", LogSpace.toString(orginalPValue), LogSpace.toString(expPValue), LogSpace.toString(newPValue));
					if (newPValue>Math.log(pvalueThreshold)) {
						Log.println(Log.Level.STANDARD, "PValue above threshold --> skipping");
						continue;
					}
					tokens.set(1, LogSpace.toString(newPValue));
				} else {
					double newPValue = orginalPValue/expPValue;
					Log.printf(Log.Level.STANDARD, "former pvalue: %e, exp. pvalue: %e --> pvalue: %e%n", orginalPValue, expPValue, newPValue);
					if (newPValue>pvalueThreshold) {
						Log.println(Log.Level.STANDARD, "PValue above threshold --> skipping");
						continue;
					}
					tokens.set(1, Double.toString(newPValue));
				}
				StringBuilder sb = new StringBuilder();
				for (String s : tokens) {
					sb.append(s);
					sb.append(' ');
				}
				Log.println(Log.Level.STANDARD, sb.toString());
			}
		} catch (IOException e) {
			Log.errorln("I/O failure while reading input file, sorry!");
			System.exit(1);
		}
		return 0;
	}	
}
